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Abstract 

This contribution presents a simple Finite Element model aimed at efficient sim- 
ulation of layered glass units. The adopted approach is based on considering 
independent kinematics of each layer, tied together via Lagrange multipliers. Val- 
idation and verification of the resulting model against independent data demon- 
strate its accuracy, showing its potential for generalization towards more complex 
problems. 
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1. Introduction 

The most frequently used transparent material in the building envelopes is 
glass. It is a fragile material, which fails in a brittle manner. This is the reason for 
using safety glasses in a situation when there is a possibility of human impact or 
where the glass could fall if shattered. 

Laminated glass is a multi-layer material produced by bonding two or more 
layers of glass together with a plastic interlay er, typically made of polyvinyl bu- 
tyral (PVB). The interlayer keeps the layers of glass bonded even when broken, 
and its high strength prevents the glass from breaking up into large sharp pieces. 
This produces a characteristic "spider web" cracking pattern when the impact 
is not powerful enough to completely pierce the glass. Multiple laminae and 
thicker glass decrease stress level, thereby increasing the load-carrying capacity 
of a structural member, too. 

The focus of this study is on the establishing a simple and versatile frame- 
work for the analysis of mechanical behavior of laminated glass units. To keep 
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the discussion compact, we restrict our attention to the linearly elastic response 
of layered glass beams in the small strain regime. The rest of the paper is orga- 
nized as follows. Methods of analysis of laminated glass beams are introduced in 
Section 2, together with a brief characterization of the proposed numerical model. 
The principles of the method are described in detail in Sections 3 and 4. In partic- 
ular, the mechanical formulation of the model is shown in Section 3. In the next 
section, the Finite Element discretization is presented. In Section 5, the proposed 
numerical technique is verified and validated against a reference analytical solu- 
tion and publicly available experimental data. Finally, Section 6 concludes the 
paper and discusses future extensions of the method. 

2. Brief overview of available methods 

The most frequent approach to the analysis of glass structural elements was, 
for a long time, based on empirical knowledge. Such relations are sufficient for 
the design of traditional windows glasses. In modern architecture there has been 
a steadily growing demand on the use of transparent materials for large external 
walls and roof systems in the recent decades. Therefore, the detailed analysis of 
layered glass units is becoming increasingly important in order to ensure a reliable 
and efficient design. 

In general, the complex behavior of laminated glass can be considered as an 
intermediate state of two limiting cases [1]. In the first case, the structure is treated 
as an assembly of two independent glass plates without any interlayer (the lower 
bound on stiffness and strength of a member), while in the second case, corre- 
sponding to the upper estimate of strength and stiffness, the glass unit is modeled 
as a monolithic glass (one glass plate with thickness equal to the total thickness of 
the glass plates). Both elementary cases, however, fail to correctly capture com- 
plex interaction among individual layers, leading to non-optimal layer thickness 
designs. Therefore, several alternative approaches to the analysis of layered glass 
structures have been proposed in the literature. These methods can be categorized 
into three basic groups: 

• methods calibrated with respect to experimental measurements [2], 

• analytical approaches [3, 4, 5], 

• numerical models typically based on detailed Finite Element simulations [6, 
7]. 

Applicability of analytical approaches to practical (usually large-scale) struc- 
tures is far from being straightforward. In particular, the closed-form solution of 
the resulting equations is possible only for very specific boundary conditions and 
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therefore have to be solved by an appropriate numerical method. Moreover, the 
analytical approaches are rather difficult to be generalized to beams with multiple 
layers. Therefore, it appears to be advantageous to directly formulate the prob- 
lem in the discretized form, typically based on the Finite Element Method (FEM). 
Nevertheless, we would like to avoid fully resolved two- or three-dimensional 
simulations, cf. [6, 7], which lead to unnecessarily expensive calculations. 

In this paper, we propose a simple FEM model inspired by a specific class 
of refined plate theories [8, 9, 10]. In this framework, each layer is treated as 
the Timoshenko beam with independent kinematics. Interaction between indi- 
vidual layers is captured by the Lagrange multipliers (with a physical meaning 
of shear stresses), which result from the conditions of inter-layer displacements 
compatibility. Such a refined approach circumvents the limitation of similar mod- 
els available in typical commercial FEM systems, which employ a single set of 
kinematic variables and average the mechanical response through the thickness 
of the beam, e.g. [11]. Unlike the proposed approach, the averaging operation is 
too coarse to correctly represent the inter-layer interactions, see Section 5 for a 
concrete example. 

3. Mechanical model of laminated beams 

As illustrated in Figure 1, laminated glasses consist mostly of three layers. A 
local coordinate system is attached to each layer to allow for an efficient treatment 
of independent kinematics. In the following text, a quantity a expressed in the 
local coordinate system associated with the j-th layer is denoted as a (,) , whereas a 
variable without an index represents a globally defined quantity, cf. Figure 1. 

Each layer is modeled using the Timoshenko beam theory supplemented with 
membrane effects. Hence, the following kinematic assumptions are adopted 

• the cross sections remain planar but not necessarily perpendicular to the 
deformed beam axis, 

• vertical displacement does not vary along the height of the beam (when 
compared to the magnitude of the displacement). 

Under these assumptions, the non-zero displacement components in each layer 
are parametrized as: 

u {i \x,z (i) ) = u ii) (x,0) + (f <i) (x)z ii \ 
w (i \x,z (i) ) = w(x), 

where i = 1, 2, 3 and z {,) is measured in the local coordinate system from the mid- 
dle plane of the z-th layer. The inter-layer interaction is ensured via the continuity 
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Figure 1 : Kinematics of laminated beam 

conditions specified on interfaces between layers in the form (i =1,2) 

1,(0 z,(!+i) 

» ( <\*,^)- M (i+ %,-^) = o. 



(1) 



The strain field in the i-th layer follows from the strain-displacement relations [12, 
11] 
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which, when combined with the constitutive equations of each layer expressed in 
terms of Young's modulus E and the shear modulus G: 



o-f(x,z w ) = E {l, ^>(x,z {l> ) and t^(jc) = G w f»{x), 
yield the expressions for the internal forces as 

Nf(x) = E»A®-^(jc,0), 
dx 



Vf(x) = kG w A (l) \<p w (x) + ^(x)\, 



dw 



dx 



M®(x) 



dx 



where b and /z (,) are the width and height of the z'-th layer, recall Figure 1, and 
k = I, A (i) = bh® a 7 W = j^bQi®) 3 stand for the shear correction factor, the 
cross-section area and the moment of inertia of the z'-th layer, respectively. 

To proceed, consider the weak form of the governing equations, written for the 
z'-th layer (the subscripts * x and » z related to internal forces and kinematics-related 
quantities are omitted in the sequel for the sake of brevity) 

L L 

f ^- (6u {i \x)) E {i) A (i) ^- (« (0 (x)) dx= f 5u (i \x)f«\x) dx + [Su(x)N (i \x)f , 
J dx dx j 

o o 

L L 

J dx kG-toA^ix) dx = J 6w(x)f®(x) dx + [Sw(x)V (i \x)f o , 
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y®(x)-<^(x)--f (w(x)) 
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to be satisfied for arbitrary admissible test fields 5u^ l \ 5(p® and 6w. In particular, 
the first three equations correspond to equilibrium conditions written for normal 
and shear forces and bending moments, respectively. The last identity enforces 
the geometrical relation (2) in the integral form, thereby allowing to treat the shear 
strain as an independent field in the discretization procedure discussed next. Fur- 
ther note that the continuity conditions (1) will be introduced directly into the 
discretized formulation, as explained in the following Section. 

4. Finite element discretization 

To keep the discretization procedure transparent, it is assumed that each layer 
of the laminated beam is divided into identical number of elements, leading to the 
discretization scheme illustrated in Figure 2. 

Following the standard conforming Finite Element machinery, e.g. [12, 11], 
we express the searched and test displacement fields at the element level in the 
form 

uf(x) * N«(x)r2„ 5uf(x) « N® (x)Sr% , 
w e (x) « N eiVV (x)r ew , 8w e {x) « N e>w (x)8r eyW , 
<pf{x) * N» (x)r® , 6<pf(x) * N» (x)6x% , 
yf(x) * H% (x)r« , Syf(x) * H%{x)5r%, 
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Figure 2: Finite element discretization of the i-th layer 



where e is used to denote the element number, * e and 8» e denote a relevant 
searched and test field restricted to the e-th element, N®. is the associated matrix 
of basis functions and r®. the matrix of nodal unknowns. In the actual implemen- 
tation, the fields u ( '\ w e and ip®, as well as the corresponding test quantities, are 
assumed to be piecewise linear. To obtain a locking-free element, the shear strain 
y® is taken as constant and is eliminated using the static condensation, see [12, 1 1] 
for additional details. 

To simplify the further treatment, we consider the following partitioning of the 
stiffness matrix K and the right hand side matrix R related to the e-th element and 
the z-th layer after the static condensation: 
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where K® =(K®) T and 

r® = [u% , u%, <p%, <p%] T , r e , w = w e , b ] T . 

Considering all three layers in Figure 2 together gives the resulting system of 
governing equations in the form 
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where the matrix A stores the nodal values of the Lagrange multipliers, associated 
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with the compatibility constraint (1), and the matrix 
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implements the tying conditions. 



5. Verification and validation of numerical model 
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Figure 3: Three point bending setup for simply supported beam 



To verify and validate the performance of the present approach, the previously 
described FEM model was implemented using MATLAB® system and compared 
against predictions of an analytical model and experimental data for a three-point 
bending test on a simply supported laminated glass beam presented in [5], see also 
Figure 3. The width of the beam is b = 0. 1 m and material data of individual com- 
ponents of the structure are available in Table 1. The glass modulus of elasticity 
is slightly lower than the conventional values of 70-73 GPa reported in the litera- 
ture and is specific for the material employed in the experiment. Moreover, as the 
PVB layer shows viscoelastic and temperature-dependent behavior, the modulus 
of elasticity corresponds to an effective secant value related to load duration of 
60 s and temperature of 22° C. 



Table 1 : Material data 





Glass PVB layer 


Young's modulus, E 
Poisson's ratio, v 


64.5 GPa 1.287 MPa 
0.23 0.4 
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Table 2 summarizes values of the mid-span deflection for a representative load 
level determined by FE-based discretization using 60 elements (30 when sym- 
metry of the problem is exploited) and the corresponding experimental values. 
Note that the discretization is sufficient to achieve three-digit accuracy in the mid- 
span deflection. In addition to the results obtained by an analytical method pro- 
posed by Asik and Tezcan in [5], the results of the analysis using ADINA® sys- 
tem and the elementary lower and upper bounds are included. In particular, the 
ADINA® model is based on the classical laminate theory, cf. [11], whereas the 
two simplified approaches assume zero or infinite compliance of the interlayer, re- 
call also discussion in Section 2. In the following discussion, e.g. 77""™ denotes the 
relative error between the numerical prediction and reference experimental value, 
while e.g. r\ an is used for the error of analytical solution when compared to candi- 
date approaches. Clearly, the results of the last three methods differ substantially 
from experimental data as well as the analytical results. The proposed numeri- 
cal model, on the other hand, shows a response almost identical to the analytical 
method, which deviates from experimental measurement by less then 6%. Such 
accuracy can be considered as sufficient from the practical point of view. 



Table 2: Comparison of results for a simply supported beam (load 50 N) 



Model 


Central deflection [mm] 




T}an 


Laminated glass beam: thickness [mm] 5/0 


.38/5 (glass/PVB/glass) 






Experiment 


1.27 




-5.2% 


Analytical model 


1.34 


5.5% 




Numerical model 


1.34 


5.5% 


0.0% 


ADINA ®(Multi-layered shell) 


0.89 


-30.2% 


-33.8% 


Monolithic glass beam: thickness [mm] 10 


(glass+glass) 






Analytical model 


0.99 


-21.8% 


-25.9% 


Two independent glass beams: thickness [mm] 5/5 (without any interlayer) 




Analytical model 


3.97 


212.6% 


196.2% 



To further confirm predictive capacities of the proposed numerical scheme, a 
response corresponding to a proportionally increasing load was investigated. The 
results appear in Tables 3 and 4. Again, the method seems to be sufficiently accu- 
rate in the investigated range of loads when considering the values of deflections 
as well as values of local stresses and strains. 



6. Conclusions 

As shown by the presented results, the proposed numerical method is well- 
suited for the modeling of laminated glass beams, mainly because of its low com- 
putational cost and accurate representation of the structural member behavior. 
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Table 3: Comparison of deflections for a simply supported beam 



Load [N] 






Central deflection [mm] 






Wexp 


Wan 


r&p [%] 


W n um 


lexp L /o J 


C m [%] 


50 


1.27 


1.34 


5.51 


1.34 


5.51 


0.00 


100 


2.55 


2.69 


5.49 


2.68 


5.10 


-0.37 


150 


4.12 


4.03 


-2.18 


4.02 


-2.43 


-0.25 


200 


5.57 


5.38 


-3.41 


5.36 


-3.77 


-0.37 



Table 4: Comparison of stresses and strains for a simply supported beam 



Load [N] 


Maximum strain [xlO 6 ] 


Maximum stress [MPa] 




^an €num 


C m [%] 


(Tan (Tnum C* [%] 


50 


112 114 


1.79 


7.23 7.34 1.52 


100 


224 228 


1.79 


14.45 14.68 1.59 


150 


336 341 


1.49 


21.68 22.02 1.57 


200 


448 455 


1.56 


28.9 29.36 1.59 



Future improvements of the model will consider large deflections and the time- 
dependent response of the interlayer and will be reported separately. 
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